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Abstract 

A technique for measuring the mean impulse response function of stationary homogeneous 
isotropic turbulence is proposed. Such measurement is carried out here on the basis of Direct 
Numerical Simulation (DNS). A zero- mean white-noise volume forcing is used to probe the tur- 
bulent flow, and the response function is obtained by accumulating the space-time correlation 
between the white forcing and the velocity field. This technique to measure the turbulent response 
in a DNS numerical experiment is a new research tool in that field of spectral closures where the 
linear response concept is invoked either by resorting to renormalized perturbations theories or by 
introducing the well-known Fluctuation-Dissipation Relation (FDR). Though the results obtained 
in the present work are limited to relatively low values of the Reynolds number, a preliminary 
analysis is possible. Both the characteristic form and the time scaling properties of the response 
function are investigated in the universal subrange of dissipative wavenumbers; a comparison with 
the response approximation given by the FDR is proposed through the independent DNS measure- 
ment of the correlation function. Very good agreement is found between the measured response 
and Kraichnan's description of random energy-range advection effects. 
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I. INTRODUCTION 



The concept of impulse response tensor of an isotropic turbulent flow lies at the heart of 
the Direct Interaction Approximation (DIA) theory, developed 50 years ago [l| by the great 
theoretical physicist Robert Kraichnan, to tackle the turbulence closure problem analytically. 
Since then, within the renormalized perturbations approach, several closure strategies have 
been proposed (a significant example is the Local Energy Transfer (LET) theory introduced 
by McComb [2|), eventually adopting a Lagrangian viewpoint, as done by Kraichnan himself 
3, 4] and others (5, 6]. In all such theories, either Eulerian or Lagrangian, closure is achieved 
by means of a closed set of integro-differential equations, where the unknowns are the two- 
points, two-times velocity correlation tensor and the response tensor itself. An exception 
is LET, where the response tensor is replaced by a renormalized propagator tensor which 
connects the velocity correlations at different times, in close analogy with the well known 
fluctuation-dissipation relation of the classical statistical physics. Recently, McComb and 
Kiyani {t| have shown how a renormalized response tensor relating the two-point covariance 
at different times can be derived; the corresponding relationship reduces to a FDR form, still 
within the theoretical framework of second-order renormalized perturbations, by introducing 
the so called time ordering approach to reconcile the time-symmetry of the correlation with 
the causality of the response. 

During the last decades, several statistics of homogeneous isotropic turbulence (HIT), ei- 
ther computed with well resolved direct numerical simulations or obtained from experiments, 
lave been compared to the corresponding theoretical predictions at increasing values of Re\ 
8], in the statistically stationary as well as in the freely decaying regime. Encouragin g re - 
sults both for the LET theory and various Lagrangian closures have been reported [q. l9l-lll|. 
Up to the present day, however, such a comparison for the impulse response function has 
never been addressed, owing to the lack of (experimental or numerical) information about 
it. Missing such a comparison is not a minor issue for Eulerian closure theories: as stressed 
in Ref. [9|, the differences among the various theoretical approaches have their roots in the 
form of the response or propagator equation, whereas the covariance equation is most often 
treated in equivalent ways. Furthermore, if the response and the two-point covariance were 
available, the degree of approximation involved in using the FDR in the context of HIT 
could be straightforwardly evaluated, indirectly gathering information about the invariant 



2 



probability distribution of the turbulent 
In recent years Luchini et al. in Ref. 



system 



12|. 



13] have proposed an original method to carry out 



an Eulerian DNS-based measurement of the mean impulse response of a turbulent flow, and 
have described the response function of a fully developed turbulent channel flow to small- 
amplitude perturbations applied at the wall. That study was conceived in the framework of 
turbulence control (hence the emphasis on wall flows and wall forcing); due to lack of isotropy, 
the response tensor is quite complicated, and does not directly relate to the previous isotropic 
theories. However, the proposed measurement technique provides us with the required tools 
to obtain the impulse response tensor for HIT, where the response function shall be intended 
to describe the response of turbulence to volume forcing. The present paper therefore aims 
at measuring the Eulerian HIT response, presenting preliminary results, obtained at low 
values of Re\, that will enable us to analyze the characteristic form and time scales of the 
response, and to compare them with theoretical predictions and assumptions. 

The paper is organized as follows. In the next fJTTJ, the definition of the impulse response 
is briefly reviewed to introduce the measurement technique, that is numerically validated 
against the available analytical viscous solution, and to discuss accuracy issues. In §1111 
the actual response function is presented and analyzed with reference to the theoretical 
background of renormalized perturbations and FDR. Lastly, §IVI is devoted to a concluding 
discussion. 



II. MEASURING THE RESPONSE FUNCTION BY DNS 
A. The definition of the impulse response function 

Following Ref. [jjj], the most general definition in wave vector space n of the instan- 
taneous impulse response tensor of a turbulent velocity field u(n, t) to an external volume 
force f(n,t), is given by the following input-output relationship between infinitesimal per- 
turbations, 6 (note the different notation from the Dirac delta function $(•)): 

bUi(K,t) = / / H in (K,K,t,t')8f n (K,t')dt'dK. (1) 

J J — oo 

It is important to underline that perturbations here assume a stochastic meaning, since 
they are superimposed to a particular random realization of u, which itself is solution of the 
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fully non-linear Navier-Stokes Equations (NSE) in Fourier space. Therefore Hi n (iz, n',t,t') 
possesses a random nature, and an integral formulation not only in time but also in wave- 
vector space is required. In fact the instantaneous response tensor plays the role of a tangent 
Green's function, related to a random and nonlinear state, and satisfies the instantaneous 
response equation: 



d_ 

dt 



+ z/k 2 ^ H in (K, K,',t,t') = 2M ijm (n) J Uj(p,t)H mn (K-p, k' ,t,t')dp+P in (K')6(K-K')8(t-t') 

(2) 

which can be derived through a stochastic Green function formalism applied to the linearized 
form of Fourier transformed NSE. In Eq. (J2D My m (/«) is the inertial transfer operator given 
by: 

M ijm (K) = -i/2(K m P ij (K) + KjP im (K)), (3) 

and Pij(n) is the projection tensor in wave- vector space, expressed as: 

Pij(K) = Sij - K~ 2 KiKj. (4) 

The locality of the response tensor in wave- vector space follows only after averaging: 

(H in )=Hin(K,t,t')5(K-K'). (5) 

Lastly, exploiting statistical isotropy and stationarity results in scalar response functions, 
respectively Q and Q, defined as follows: 

U in {K, t, t') = P in {n)g{K, t, 0, (6) 

G(k,t) = G(K,t,t-T). (7) 
The causality property holds for both the previous functions, hence : 

Q(k, t) = for r < and V/c. (8) 

This is obviously a consequence of the realizability of the dynamical system that is being 
described through its impulse response. As indicated by Kraichnan the scalar response 
is a real, unit-bounded function: 

\Q(k,t)\ < Q(k,0 + ) = 1, VrXDandV/t. (9) 
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FIG. 1: Compensated energy spectrum for HIT: the function E(k) computed with the present 
DNS code at several values of Re\ is compared with results from Ref. [3] at Re\ = 84. 

B. The Direct Numerical Simulation 



The measurement of Q described in this paper is carried out by means of a forced DNS 
of stationary HIT on a cubic domain, whose edge length L is chosen to be L = 2tt for 
convenience, so that the fundamental wave number is k = 2n/L = 1 without loss of 
generality. A numerical code has been developed on purpose and equipped with parallel 
(shared-memory) computing capabilities. The code implements a classical Galerkin-Fourier 
scheme applied to the velocity-vorticity formulation of the incompressible Navier-Stokes 
equations. In the present context, this formulation presents interesting advantages in terms 
of memory requirements. Exact removal of the aliasing error is obtained with the 3/2 
zero-padding rule; time integration is carried out by means of a third-order low-storage 
Runge-Kutta (Williamson) scheme; see Ref. HQ for additional numerical details. The 
forcing scheme has been carefully implemented following the provisions stated in Ref. 15], 
from which the notation adopted below is borrowed. The Kolmogorov scale is indicated with 
T], with = t]^ 1 , the instantaneous dissipation rate is e, the forcing- containing shell is Kf 
and the mean energy injection rate is P, that equals (e) at statistical stationarity. Then the 
adopted feedback-acceleration forcing |l5j is formulated in wave number space as follows: 



t^-^m 1 ^ < 10) 
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where kf(t) represents the kinetic energy of the modes within the forced shell Kf and h(K; Kf) 
is the related indicator function: 

{1, \K\ < K f , 
(11) 
0, otherwise. 

A standard resolution of K max r] = 1.5 is adopted, where K max indicates the maximum re- 
solved wave-number in each direction of the Fourier space. The numerical code has been 
thouroughly verified by running conventional simulations of stationary HIT. The computed 
energy spectra at various Re\ compare very well to available results. A comparison of this 
kind is shown in Fig. [TJ that shows excellent agreement between our computed energy spec- 



tra and those pubblished in Ref. [15j. The spectral code has been run on a machine equipped 
with 4 Opteron 2378 processors, where a case with N = 256 has a memory requirement of 
940MB and a typical execution time of 11 seconds for one Runge-Kutta time step. 

C. The response measurement technique 
In Ref. 



13j Luchini, Quadrio & Zuccher propose an innovative method for measuring the 
linear impulse response of a turbulent velocity field, resorting to the statistical statement of 
the input-output relation for a linear system, i.e. the input-output correlation. This approach 
is primarily motivated by the problem of low signal-to-noise ratio (S/N) that one would 
face, should the response function be measured according to its definition. Indeed the linear 
response of a non-linear dynamical system is obtained by means of infinitesimal perturbations 
around an equilibrium state, which has a stochastic meaning in the description of turbulence. 
Therefore impulsive perturbations externally introduced into the turbulent field to measure 
its linear response must be extremely small compared to the natural turbulent fluctuations 
for Eq. (pQ) to hold; as a consequence, their effect is buried into turbulent noise. 

By definition the impulse response is the output of a linear system when either harmonic 
or impulsive signals are used as inputs. However, for a linearized turbulent system, the use 
of a proper statistical probe instead of a deterministic one will dramatically improve the 
computational efficiency of the overall measurement procedure. This is the case of using a 



white-noise process in input to the system. Indeed it is well known from filtering theory 18] 



that when a linear system is fed with white noise, the correlation between the input and the 
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output is proportional to the impulse response of the system, owing to the delta-correlated 
property of the white-noise process. We employ an externally generated random volume 
forcing as the input; by computing its cross-correlation with the velocity field, the whole 
wave-number dependency of the response function is obtained at once. At the same time, 
forcing is uniformly distributed over time and space, thus leading to improved S/N and 
larger allowed amplitudes within the linearity constraint. Therefore this strategy performs 
much better than a deterministic forcing, be it either harmonic or impulsive, that would 



lead to computationally unaffordable simulations, as highlighted in Ref. 13 1. 
Starting from Eq. (jTJ), the input-output correlation can be written as: 



p p+oo 

(&Ui(K,t)6f j (-K,t-T))= / / H ln (K,t-t')5(K f -K) (SfniK^t^bfji-^t - T)} dt'dK,', 

J J — oo 

(12) 

where Eq. (J5J) has been used owing to the average operator, and the response causality 
property allows the extension towards +00 of the upper bound of time integral. Assuming 
5fj(n,t) = ewj(n,t), being e e IR + a scale factor and Wj(K,t) an independently generated 
zero-mean white-noise field with identity covariance matrix: 

(6/ n («, 1?)6fj(-K, t-r)) = e 2 5 nj 5(t' -t + r), (13) 

the cross-correlation at the l.h.s. of Eq. ffT2]) will result in the properly scaled response 
tensor: 

(SuifotjS/^-M-r)) = e 2 ^(K,T). (14) 

We shall denote by u(k, t) the turbulent velocity field when volume forcing with white 
spectrum is applied. If the perturbation is small enough for linearity to hold, i.e. e 1, it 
follows that: 

u(k, t) = v(k, t) + 5u(k, t), (15) 

where v(n,t) indicates a different realization of the turbulent fluctuating field respect to 
the original field u(n,t), as a consequence of non-linearity and stochastic behavior of NSE. 
Then computing the correlation between u and 5/ results in: 



^1 — T)) = \ [(vi(K, tM(-K, t-r))+ (M*. t - r)>] . (16) 
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Since the applied random perturbation on forcing is uncorrelated to turbulent fluctua- 
tions, the term (vi(n, t)bfj(—K, t — r)) will be averaged out in the previous equation, leading 
to: 

(u i (K,t)5f j (-K,t-r)) 



€ 2 



Hufcr), (17) 



where the input-output correlation law, Eq. (1121) . has been used to handle the non- vanishing 
term (the second term) at r.h.s. of Eq. ([To 7 ]) . In this way it is still possible to measure the 
turbulent response using the cross-correlation between the white-noise input and the whole 
turbulent velocity field. 

At this point it may be useful to note explicitely that no relation exists between the 
energy- driving forcing, Eq. ( TTUj) . and the white- noise forcing applied for the response mea- 
surement, Eq. ( fi~3l) . The former obviuosly represent a mere artificial, but unavoidable, 
mechanism to drive the flow and maintain statistical stationarity, via supply of energy at 
large scales. The white noise, on the other hand, is an external field of volume force that 
enters the stationary turbulent system, which includes the energy-driving forcing. The j-th 
component of the white noise field, Wj(K,t) is defined as: 

Wj(n, t) = exp(z27T0), (18) 

where <j) is the output of a random number generator with an uniform probability distribution 



in the interval [0,1] 32|. Hence the white noise is independent from both the turbulent 
fluctuations and the energy- driving forcing applied to them through the feedback formula 
(JTDJ). A feedback forcing loop should be considered when looking at the linear response 
for wave numbers contained in the forced shell, but, as already discussed, these are not of 
physical interest. 

In the HIT case Eq. ([6]) provides us with a convenient way of accumulating just the 
simple scalar version of the response tensor, by means of shell averaging over tensor trace: 

Hu(k, t)cIS(k) = Stxk 2 Q(k, t). (19) 

When measuring Q(k, t), a proper spatial and temporal discretization must be adopted. 
While in a DNS the discretizazion in k is easily derived from the Fourier representation of the 
velocity field (as for the energy spectrum E(k), see Ref. [jjj]), the definition of the r-step is 
less obvious. Both the time resolution of the white-noise delta correlation, At w , i.e. the time 
interval between successive updates of the random numbers, and the averaging time T av , i.e. 
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the time interval over which statistics are computed, must be chosen such that Q is properly 
described and at the same time the computational requirements of the numerical simulation 
are kept reasonable. If r m j ra indicates the smallest time scale at which proper convergence of 
the response is sought, At w must be chosen so that At w < r min . Assuming uniform sampling 
of the response in N c time instants separated by At, the time horizon available to represent 
the decay of the entire response must be greater than the whole response decay time r max at 
the lower wave number in the range of interest, i.e. N c Ar > r max . Proper convergence of the 
average response obviously requires T av /T max 3> 1. Indeed, while Kd controls At resolution 
and then At, i.e. the time integration step, the largest inertial wave number dictates N c and 
T 

Given such contrasting requirements, characterizing the function Q(k,t) in the whole 
universal range of scales via a sole measurement is possible but computationally demanding. 
Hence the entire function Q(k,t) can be measured through more than one uniform r-grid, 
so that the response is probed within several sub-ranges of scales, leveraging their reduced 
extent. Q(k,t) is then measured in a wide range of scales via a limited number of DNS 
runs, each of which requires roughly the same computational effort. These simulations 
are independent, and can be run simultaneously if the available computing power allows. 
However, for the results to obey the linearity costraint, the level of the introduced "noise 
energy", eAr w , must be kept constant across the different At resolutions adopted at different 
scales. This means that a larger At implies a reduced noise amplitude e and a longer 
averaging time. This is partially compensated by the larger time step size allowed by the 
time resolution of the response at lower wave numbers. 

D. A test case: the purely viscous Stokes' response 

The Stokes or viscous response represents the zero-order term in the expansion series of 
Q as introduced in the context of renormalized perturbations, see Ref. [9I. Tioj] . The Stokes 
response, 

£(0) 

, can be easily derived from Eq. ([2]) after removal of the non-linear terms, thus 
providing the solution for pure viscous dynamics of the velocity field. Its analytical form 
reads: 



It is important to notice that the Stokes response has a deterministic nature, owing to 




(20) 
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TABLE I: Parameters for the DNS of HIT carried out in the present work. 



N 


l^max/ ^0 




P 


128 


42 


28 


1 


192 


63 


42 


1 


256 


84 


56 


1 



Kf/Ko Re=\^—j Re\ u y — 

3 20 55 1.7862 

3 34 77 1.8453 

3 49.5 94 1.8611 



the linearity of the Stokes operator: Kraichnan usually refers to it as "statistically sharp". 
The exact Stokes solution provides an useful tool for the validation of the full measurement 
procedure. To this purpose, the Stokes response can be retrieved from a DNS of the fully 
non-linear NSE through a numerical linearization. In this way the algorithm employed for 
the measurement in the turbulent case is exactly that previously described in §11 CI but a 
null initial condition is adopted, the energy-driving forcing of Eq. (fTUj) is turned off, and only 
the white-noise perturbation is applied. If e 1, no evolution toward turbulence dynamics 
is produced, and non-linear terms C(e 2 ) can be neglected with respect to the linear ones 
(9(e) defining the Stokes equation. 

The Stokes response has been measured in numerical experiments with a spatial resolution 
of 32 3 modes (before dealiasing) and N c = 50. In these simulations At = Ar w is adopted, 
and several values of At w are used to investigate time resolution effects. In Fig. |2] (top and 
center) the time decay of the Stokes response at k/k q = 8 is plotted. The exact solution and 
the measured one agree very well (top figure). At small viscous time separations, tuk, 2 < 1, 
proper converge of the measured response towards the exact one is observed (center figure) 
to depend on the At w resolution. Convergence of the response for different noise amplitude 
e and averaging time T av has been additionally verified (not shown). Lastly, the measured 
Stokes response plotted at different wave numbers (bottom figure) is observed to possess the 
expected collapse when local viscous time scaling is adopted for r. 



III. THE TURBULENT RESPONSE 



Several DNS have been run to measure the impulse response of homogeneous isotropic 
turbulence. They are grouped into simulations with 128 3 , 192 3 and 256 3 Fourier modes, be- 
fore dealiasing. Table [J summarizes the discretization parameters of the simulations carried 
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FIG. 2: Time decay of the measured Stokes response £ (0) at different At with e = 0.001 and 
T av vK 2 = 733.76. Top: comparison between the measured and exact at fixed k/kq = 8. 
Center: zoom of top plot for tvk? <C 1. Bottom: Stokes response at several wave number k/kq vs. 
non-dimensional time separation, emphasizing collapse with local viscous time scaling run 2 . 

out without white-noise forcing, whereas Table [III lists all the simulations run to measure 
the response function, together with the values of the parameters used to discretize and 
measure Q(k,,t). The attained values of Re\ are low or moderate, ranging from Re\ = 55 
to Re\ = 94. Moreover, given our limited computational resources, the response is probed 
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TABLE II: Discretization parameters for the DNS-based measurement of the response function. 



N Re x 


Run 


N c 


ArUnKri 


e 


r\_u U (l 




1 


150 


0.05202 


0.00093 


1.3004e5 


128 55 














2 


150 


0.0322 


0.0015 


5643.7 




1 


150 


0.0484 


6.667e-4 


9680 


192 77 














2 


150 


0.0322 


0.001 


5643.7 




1 


150 


061 A 


4.3529e-4 


9981.6 


256 94 


2 


150 


0.0376 


7.1154e-4 


6580 




3 


150 


0.0267 


0.001 


2069.2 


only at high wave- 


-numbers. As explained in 


511 C[ lone averaging 


times are in fact required 


for proper convert 


fence of the 


mean reponse 


at low wavenumbers. 


However, 


it is important 



to emphasize here that the proposed measurement technique is capable of measuring the 
impulse response at any scale, provided that adequate computational resources are available. 

In §111 Al the linear behavior and the convergence of time averages are demonstrated, 
whereas in §111 Bl the behavior of the response and its scaling within the dissipative universal 
subrange are investigated. The correlation function is also computed during the stationary 
reference HIT simulations, so that a comparison with the impulse response function in terms 
of the classical FDR will be given at least in the range of scales here considered. 

A. Linearity and time average 

A key issue when measuring the response function Q{k, t) is the proper choice of ampli- 
tude e for the white-noise forcing. Indeed the true turbulent impulse response reduces to 
its linear counterpart Q only for vanishing "noise energy", i.e. when eAr w — > 0. In a finite 
setting, a reasonable preliminary requirement is that the white-noise forcing does not affect 
turbulence statistics appreciably. Then, suitable convergence of the measured response to 
Q is observed when the function Q/e becomes indipendent on e. Indeed, as discussed in 
§11 C\ given a time resolution At w , e represents the linearity control parameter. Since the 
white-noise forcing is spatially distributed over all the scales, the linearity threshold is fixed 
by pertubations effects on smallest scale dynamics, i.e. the viscous scales: when looking for 
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1.4 



FIG. 3: Comparison between compensated energy spectra from HIT and from white-noise- forced 
HIT employed for measuring G(k,t), Re x = 94 {N = 256) (cfr. Table HJ). Note the absence of 
white- noise effects for k/k^ < 1. 

the response at lower wave numbers with larger At = At w , e must be reduced to preserve 
the linear response at higher wave numbers. 

Fig. [3] provides a comparison between the energy spectrum E(k) computed in standard 
HIT DNS and those from simulations forced with white noise at Re\ = 94. Analogous results 
(not shown here) holds for the other spatial resolutions listed in Table HT1 The value of e, that 
should be maximized in order to increase S/N and hence to reduce the required averaging 
time, is chosen so that marginal effects on the spectrum are confined within the numerical 
wavenumbers larger than the Kolmogorov scale, k > k^. Moreover, table II III quantifies the 
little variations in statistics like (e) and (k) due to white-noise forcing: A (e) and A (k) are 
less then 0.4% and 1.9% respectively. A(e), that is computed with respect to the exact 
asymptotic value P, is of the same order of the variations of (e) observed in different runs of 
standard HIT DNS. A (k), that is computed against the value of (k) obtained from standard 
HIT DNS (Run 0), seems to be relatively larger. However, it should be recalled how the 
accurate convergence of this statistic, that belongs to large scales, requires averaging over 
many turnover times. Indeed, the observed A (k) are of the same order of (k) fluctuations in 
standard HIT DNS, under the feedback action of the energy-driving forcing scheme with an 
averaging time of T av (PK 2 f y/ 3 w 250, to which standard HIT averaged values are referred. 

The convergence of the measured Q{k, t) with respect to the averaging time T av is also 
verified. In Fig. H] responses measured with increasingly larger values of T av are shown 
for Run 2 at Re\ = 94. Adequate convergence is obtained at representative wavenumbers 
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TABLE III: Effect of the white-noise forcing on (e) and (k), for various spatial resolutions. The 



reference simulations without whit 


e noise are 


indicated as "Run 0". 






N Re x 


Run 


(e) IP 


(k)/(P/ Kf )^ 


A (e) % 


A (k) % 







0.999595 


4.7856 


0.0405 


- 


128 55 


1 


1.00057 


4.7904 


0.057 


0.1 




2 


1.00056 


4.7374 


0.056 


1.007 









o.yjuoo 


D D91 9 




192 77 


1 


1.00177 


4.9612 


0.177 


1.83 




2 


1.U025 


5.0063 


0.25 


0.94 







1.00155 


5.1958 


0.155 






1 


1.00268 


5.1235 


0.113 


1.392 


256 94 














2 


1.0007 


5.1082 


0.085 


1.686 




3 


1.00354 


5.2561 


0.354 


1.161 





T /N Ax - 6.5e+002- 

av c 




- - - T /N Ax = 7.8e+002 

av c 




T /N Ax = 9.2e+002 

av c 




T /N Ax = 1.2e+003 

av c 


\ 













FIG. 4: Convergence of the measured Q(k,t) with respect to the averaging time T av . Results for 
Run 2 at Rex = 94 (cfr. Table ITi]) are shown at rapresentative wavenumbers k/kj = 0.75. 

ft/ftd — 0.75 when T av j (N c At) > 920. At larger averaging times, the response curves become 
indistinguishable. A similar behavior has been verified for the other numerical experiments 
reported in Table [ITI 
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FIG. 5: Measured response function from Run 3 of Table ITU at Re\ = 94. Top: time decay of the 
response function at the Kolmogorov scale (n/nd = 1) 5 compared with the DIA solution Eq. (|21f) 
and the viscous Gaussian-convective solution Eq. (|22[) . Center: zoom of the top plot for tkuq <C 1. 
Bottom: zoom of the top plot at large tkuq. 



B. The response function and its scaling in the viscous universal subrange 



The response function measured via the procedure illustrated above is first compared 
with its available analytical approximations, as given in the original DIA theory, see Ref. 
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Q: 



2 ^Ji(2u kt) 



Q(k,t) = exp(—vn r) , (21) 



and in the analysis of random convection effects [9|, |20[ from which the viscous Gaussian- 
convective response Qgc{k,t) can be introduced: 

Ggc(k, t) = exp(— vk 2 t — -u^k 2 t 2 ), with r > 0. (22) 

In both the previous equations Uq represents the r.m.s. value of turbulent fluctuations and 
it can be easily recognized the Stokes term, Eq. ( 1201) . which reflects the viscous response of 
the corresponding linear operator. It is important to recall that while the non-viscous term 
of Eq. (l2Tj) is derived as an approximated solution to the DIA equations, the corrisponding 
one of Eq. fl22|) empirically follows from the analogy with the solution of the idealized 
problem of pure random convection introduced by Kraichnan in Ref. {20] with the Random 



Galilean Invariance (RGI) postu 



inertial-range scaling. Refs. 21 



ate to explain the failure of DIA in yielding a Kolmogorov 



22] give a more recent investigation on the role of random 



convection effects and RGI in renormalized perturbation expansions of the NSE. 

A comparative view of these three response functions at n/nd = 1 is provided in Fig. 
[5] for Re\ = 94, Run 3. At time separations smaller than the local energy time scale, i.e. 
for tuqk < 1, the true measured response is in good agreement with the DIA response 
function and the viscous Gaussian-convective solution. The latter result does not come as a 
surprise: even though the turbulent field is definitely non-Gaussian, at times smaller than 



the characteristic correlation time the Gaussian approximation still applies, see [19] . The 
unexpected result, however, is that the Gaussian convective solution still approximates very 
well the measured response at larger times, whereas the DIA solution clearly deviates from 
it. 

This evidence provides further motivation for investigating the convective response scaling 
in the viscous universal subrange. Response functions rescaled accordingly are plotted in 
Fig. |6j For the entire range of values of Re\ considered in the present work, the convective 
scaling of the response function is clearly assessed. 

To further support this statement, Fig. [7] shows the response functions tentatively plotted 
with Kolmogorov viscous scaling: it is evident that such scaling does not produce as good a 
collapse of the different curves when compared to the convective scaling employed in Fig. EJ 
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FIG. 6: Measured response functions plotted with convective scaling at different values of Re\. 
The responses are plotted versus non-dimensional time separation tuok, at several wavenumbers 
within the universal dissipative subrange (cfr. Table HI]) . Top: Re\ = 77. Bottom: Re\ = 94. For 
better clarity response functions are plotted using one of every two of the N c values. Convective 
scaling produces a good collapse of the curves. 

C. The correlation function and the FDR 

The mean correlation tensor is here introduced directly in its spectral form: 

Qij{K,t,t') = (Ui(K,t)Uj(-K,lf)) , (23) 

which reduces to the scalar function Q(k, t) in the homogeneous, isotropic stationary case: 

Q ij (K,t,t , ) = P ij (K)Q(K,t-t'). (24) 

The correlation function Q(n, r) has been computed thanks to the DNS simulations 
carried out without white-noise forcing. At the various values of Re\ considered, the cor- 
respondent smallest At time resolution employed for measuring the response function has 
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FIG. 7: Measured response functions with Kolmogorov scaling at Re\ = 94. The responses are 
plotted versus the non-dimensional time separation r(P^) 1 ^ 3 , at several wavenumbers within the 
universal dissipative subrange (cfr. Table HI)) . For better clarity response functions are plotted 
using one of every two of the N c values. Kolmogorov scaling is not successful in producing a 
collapse of the curves. 
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FIG. 8: Measured correlation function plotted with convective scaling. The normalized correlation 
function is plotted versus non-dimensional time separation tuqk, for several wavenumbers in the 
universal dissipative subrange at Re\ = 94. Convective scaling produces a good collapse of the 
curves. 



been used. Discretization details can be found in Table HIl The number of time separations 
at which the correlations are stored is N c = 200 for the two cases respectively at Re\ = 55 
and Re\ = 77, while N c = 175 has been employed for the case at Re\ = 94 due to memory 
limitations. Proper convergence of the results with respect to T av has been verified, according 
to what has been done for the response function itself. Similarly to what has been observed 
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FIG. 9: Measured correlation function plotted with viscous scaling. The normalized correlation 
function is plotted versus non-dimensional time separation t(Pk^) 1 ^ 3 , for several wavenumbers in 
the universal dissipative subrange at Re\ = 94. Viscous scaling does not produce a collapse of the 
curves. 

for Q(k,t), the scaling of the normalized correlation function, Q(k, t)/Q(k, 0), is captured 
by the local convective time scale (km ) _1 in the universal viscous subrange investigated 
here. This is illustrated for Re\ = 94 in Fig. [HI whereas the inadequacy of Kolmogorov 
viscous scaling is shown in Fig. |9j The same behavior can be also observed to hold for the 
correlations obtained at Re\ = 55 and Re\ = 77 (not shown here). 

The response and correlation functions, as measured from our DNS experiments, can be 
compared through the well known FDR: 



Q(k,t) = G{k,t)Q(k,0). 



(25) 



This relation has been originally derived in the context of Hamiltonian dynamical sys- 
tems at equilibrium for which a canonical distribution holds. Only in the last decades, the 
applicability of the FDR to the wider class of non-linear chaotic dynamical systems has been 

25I] . For this class of systems (to which fluid turbulence 



addressed on a theoretical basis 
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belongs) a generalized FDR is demonstrated to hold, provided the system is dynamically 
mixing: only when a Gaussian distribution holds for the invariant probability distribution, 
the generalized FDR reduces to the classical form of Eq. (I2"5l) . Obviously Eq. fl2"5]) cannot be 
exact for fully developed fluid turbulence, for which both experimental and numerical inves- 
tigations have shown marked departures from Gaussianity, with long tails in the PDF and 
intermittent behavior. However on an intuitive ground, one would expect a proportionality 
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between the response and the correlation function to hold, at least in terms of characteristic 
time scales, respectively indicated by tq{k) and tq(k) [33]. FDR has been then successfully 
applied in the context of climate study on sensitivity analysis with respect to external per- 



in Refs. 
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turbations and parameters |26l . |27| , as well as in viscosity renormalization [28| . Nevertheless, 
2 7| | it is noted how in many such attempts the Gaussian form of the FDR has 



been often acritically invoked, with little awareness about its inherent limitations. More- 
over, in the field of spectral closures, different opinions exist on the possibility to recover 
the classical FDR in the Eulerian rather than in the Lagrangian framework. In Ref. [5:] the 
Gaussian form of the FDR is exactly recovered within the Lagragian renormalized approxi- 
mation of turbulence, where the Lagrangian response function is introduced. In the Eulerian 
frame the proper use of the FDR has been recently addressed by Kiyani and McComb Q. 
In their paper they show that FDR as stated in Eq. (1251) is exact up to second order in 
renormalized perturbation expansions of NSE, hence it can be properly used in related clo- 
sure formulations {22 1. However, Kraichnan suggested 30] that even a valid Gaussian FDR 
would not immediately be a step forward in the closure problem. In Kraichnan's view, the 
strong departure from equipartition in the inertial range is not followed by a corresponding 
strong violation of the Gaussian FDR in the Eulerian frame. This is because large-scale 
random convection dominates the decay of both the response and the correlation functions, 
with corresponding time scales for mode k ruled by the local characteristic convective time 
(fcuo) -1 - Kraichanan's analysis thus implies that the local dynamics cannot be captured 
by the elementary FDR, and the expected deviations can be found only by looking to a 
generalized FDR, that involves a Lagrangian form of the statistics. 

These considerations motivate investigating the approximation introduced by the Gaus- 
sian FDR within the Eulerian frame: a preliminary assessment is given Fig. [10] where 
G(k, t) and t)/Q(k, 0) are plotted together for k fixed at the Kolmogorov scale, i.e. 
for K/Kd = 1. As expected from theoretical arguments, see for example Ref. 19|], a longer 
decorrelation time is observed for Q(n, t)/Q(k, 0) when compared to Q(k,t). However the 
time scales rg(n) and tq(k) turn out to be of the same order. At fixed Re\, the plots in 
convective units of Fig. [10] (top and center) suggest that the response and the correlation 
functions are strictly related within the whole dissipative subrange of scales, owing to their 
inherent energy-convective scaling property previously discussed. 

When examining the response and the correlation functions obtained at different values 



20 



of Re\, one first observes that, in agreement with the very good approximation provided by 
the analytical viscous Gaussian-convective formulae of Eq. fl22|) to the response function, 
the latter is well described as an universal function of the adimensional variable tkuq. The 
same observation does not hold true for the normalized correlation function Q(k, t) / Q(k, 0) 
which once plotted in convective scaling shows a residual dependence on Re\. In particular 
when increasing Re\, the correlation function moves towards the response function, and 
this implies that the approximation involved by the classical FDR, Eq. (125]) . is gradually 
improving. For completeness the two functions t) and Q(n, r)j Q(k, 0) are also plotted 
in terms of Kolmogorov viscous units, as shown in Fig. [TO] (bottom). When this scaling is 
employed, neither the response nor the correlation show a collapse. 



IV. CONCLUSIONS 



The mean linear response function of homogeneous isotropic turbulence to an impulsive 
body force has been measured through a number of numerical experiments carried out with 
DNS, at low and moderate values of the Reynolds number Re\. The measurement method 
leverages a white-noise forcing to probe the flow within the linearity constraint while man- 
taining the computational effort at reasonable levels. The method employed for measuring 
the response of the full turbulent flow has been thouroughly validated by computing the 
response function for purely viscous dynamics: the very same procedure yields this simpler 
response, for which an exact analytical expression is available to compare with. Based on this 
test case, the proper convergence with respect to the parameters of the time discretization 
has been verified. The methodology proposed here for measuring the response function has 
then proved effective in the quantitative description of the whole time decay of the response 
within the universal equilibrium range of scales. Our results have been verified both in terms 
of linearity of the response with respect to the amplitude of the forcing and adequateness of 
the time averaging. The same direct numerical simulations have been additionally employed 
for determining the turbulence correlation function. Its examination within the same range 
of scales has allowed us to preliminarly address the approximations involved by the classical 
fluctuation-dissipation relation when applied to turbulence dynamics. 

The analysis of the response function in the universal dissipative subrange confirms the 
theoretical prediction of energy-convective scaling for both the response and the normal- 
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FIG. 10: Measured response functions Q(k,t) (continuous lines) and measured normalized correla- 
tion functions Q(k, t)/Q(k, 0) (dashed lines) at the Kolmogorov scale k = for several values of 
Re\. Top: convective scaling. Center: convective scaling, but using only turbulent-diffusive part 
of the response is used. Bottom: viscous Kolmogorov scaling. 

ized correlation functions, and, as shown by Figs. [6] and [H establishes such scaling as the 
dominant one, at least in the rather limited range of Re\ considered here. A somewhat 
surprising result is that the analytical solution provided by Kraichnan in Ref. 20j to the 
problem of idealized convection turns out to be an extremely good approximation of the 
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measured response function, with small deviations limited to the tail region, as shown in 
Fig. E 

When comparing the normalized correlation function and the response function, a longer 
decorrelation time is observed for the former, as suggested by the theoretical arguments put 
forward in Ref. [19|]. Both Q(k,t) and Q(k, t)/Q(k, 0) obey the same convective temporal 
scaling within the dissipation range, hence the two time scales tq(k) and tq(k) are in an 
approximately constant ratio. Obviously the Gaussian form of the FDR, Eq. (125]) . is not 
exactly satisfied, as witnessed from the departure between r) and t)/Q(k, 0) in 
Fig. QUI Nevertheless, Eq. f l23|) remains a good approximation in terms of characteristic time 
scales, even at the moderate value of Re\ considered here and in the dissipation subrange, 
where less agreement would be expected in comparison to the inertial subrange, which 
is the proper context in which the FDR should be considered [311 ] . Moreover, the FDR 
approximation in the present range of scales appears to be increasingly better supported 
when the value of Re\ is increased, as shown in Fig. [TU1 This last conclusion is in partial 
agreement with a previous study by Biferale et al. [is] , who examined the response function 
within the inertial range of scales as extracted from the shell model. In that work the concept 
of halving-time statistics was introduced to better characterize the time properties of the 
response function for lower shells, where the proper r-convergence of the response cannot be 
easily achieved. The ratio between characteristics times is still constant in the inertial range, 
but tq(k) and tq(k) show Kolmogorov inertial time scaling. Both Kraichnan's arguments on 
random convections effects as well as the more recent and related discussion on the validity 
of the FDR in the context of turbulence, Ref. 30], are strongly supported by present results. 
However Kraichnan indicates that the dominance of energy-advection effects on both the 
Eulerian response and the correlation functions are expected to extend to the dissipation 
range only at high values of Reynolds number, while this has been found in the present 
work to happen already at low or moderate values of Re\ addressed here. One possible 
explanation might be provided by considering the energy-convection effects as a feature of 
turbulence that remains limited to the dissipative subrange of scales, so that the presence of 
significant scale separation from the energy scales would let the random convection picture 
to hold: however the same could not be true for inertial scales at higher Re\. 

A more thorough description of the response function and of its relevant time scales, 
together with a precise assessment of the approximations involved by the classical FDR, 
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obviously call for an extension of the present study towards much higher values of Re\, 
so that a well-defined inertial range can develop. When such data will be available, the 
question about a possible asymptotic vanishing of the convective scaling in favor of a true 
Kolmogorov scaling could be properly answered, thus enlightening the framework of Eulerian 
closure theories. To this purpose, an analysis using halving-time statistics can be exploited 
to accurately characterize the properties of the response function in time over a wide range 
of scales. If Kolmogorov scaling will indeed be recovered at higher Re\, then the local 
relaxation processes of the turbulent response would be captured, opening a new scenario 
in the understanding of turbulence physics and modeling. 

Acknowledgments 

The authors would like to thank Dr. F. Martinelli for suggesting the Stokes test case 
described in §11 Dl and for the helpful discussions. We gratefully acknowledge the use of the 
computing system located at the University of Salerno and the discussions with Prof. P. 
Luchini. 



[1] Kraichnan R. H., J. Fluid Mech. 5, 497 (1959). 

[2] W. McComb, J. Phys. A: Math. Nucl. Gen. 7, 632 (1974). 

[3] Kraichnan R. H., Phys. of Fluids 7, 575 (1964). 

[4] Kraichnan R.H., J. Fluid Mech. 83, 349 (1977). 

[5] Y. Kaneda, J. Fluid Mech. 107, 131 (1981). 

[6] S. Kida and S. Goto, J. Fluid Mech. 345, 307 (1997). 

[7] K. Kiyani and W. McComb, Phys. Rev. E 70, 066303 (2004). 

[8] T. Ishihara, T. Gotoh, and Y. Kaneda, Annu. Rev. Fluid Mech. (2009). 

[9] W. McComb, The Physics of Fluid Turbulence (Oxford University Press, 1990). 

[10] W. McComb, M. J. Filipiak, and V. Shanmugasundaram, J. Fluid Mech. 245, 279 (1992). 

[11] Kaneda Y., Phys. Fluids 5, 2835 (1993). 

[12] L. Biferale, I. Daumont, G. Lacorata, and A. Vulpiani, Phys. Rev. E 65, 016302 (2002). 

[13] P. Luchini, M. Quadrio, and S. Zuccher, Phys. Fluids 18, 1 (2006). 



24 



[14] P. Sagaut and C. Cambon, Homogeneous turbulence dynamics (Cambridge University Press, 
2008). 

[15] A. G. Lamorgese, D. A. Caughey, and S. B. Pope, Phys. of Fluids (2005). 

[16] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral methods - fundamentals 

in single domains (Springer, 2006). 
[17] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods. Evolution to 

Complex Geometries and Applications to Fluid Dynamics (Springer, 2007). 
[18] A. H. Jazwinski, Stochastic processes and filtering theory (Academic Press, New York, 1970). 
[19] Leslie D. C, Developments in the theory of turbulence (Clarendon Press, Oxford, 1973). 
[20] Kraichnan R. H., Phys. of Fluids 7, 1723 (1964). 

[21] W. McComb, V. Shanmugasundaram, and P. Hutchinson, J. Fluid Mech. 208, 91 (1989). 
[22] W. McComb, Phys. Rev. E 71, 037301 (2005). 

[23] M. Falcioni, S. Isola, and A. Vulpiani, Phys. Lett. A 144, 341 (1990). 

[24] M. Falcioni, S. Isola, and A. Vulpiani, Phys. Fluids A 3, 2247 (1991). 

[25] U. Marconi, A. Puglisi, L. Rondoni, and A. Vulpiani, Phys. Reports 461, 111 (2008). 

[26] T. L. Bell, J. Atmos. Sci. 37, 1700 (1980). 

[27] G. Lacorata and A. Vulpiani, Nonlin. Processes Geophys. 14, 681 (2007). 
[28] G. F. Carnevale and J. S. Frederiksen, J. Fluid Mech. 131, 289 (1983). 
[29] W. McComb and K. Kiyani, Phys. Rev. E 72, 016309 (2005). 
[30] R. H. Kraichnan, Physica A 279, 30 (2000). 

[31] Rosenblatt M. and Van Atta C, eds., Statistical models and turbulence (Springer Verlag, 1972), 

vol. 12, chap. Comparison of some approximations for isotropic turbulence, pp. 148-194. 
[32] Note that a white noise field with identity covariance matrix is used in the definition of 

5/j(ft, t), Eq. (fT3|) : the identity covariance scaling factor for Wj(n, t) could be easily embedded 

in the definition of e without losing of generality. 
[33] tq(k) and tq(k) are to be intended respectively as the separation time r by which Q(k, t) and 

Q(k, r) reduce to the same percentage of their initial value. 



25 



